{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We covered the idea of simulating an arbitrary transfer function system in [a previous notebook](../1_Dynamics/5_Complex_system_dynamics/Simulation%20of%20arbitrary%20transfer%20functions.ipynb). What happens if we have to simulate a controller (specified as a transfer function) and a system specified by differential equations together?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Nonlinear tank system\n",
    "\n",
    "Let's take the classic tank system, with a square root flow relationship on the outflow and a nonlinear valve relationship.\n",
    "\n",
    "<img src=\"../assets/tanksystem.png\">\n",
    "\n",
    "\\begin{align}\n",
    "      \\frac{dV}{dt} &= (F_{in} - F_{out}) \\\\\n",
    "      h &= \\frac{V}{A} \\\\\n",
    "      f(x) &= \\alpha^{x - 1} \\\\\n",
    "      F_{out} &= K f(x) \\sqrt{h} \\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = 2\n",
    "alpha = 20\n",
    "K = 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initial conditions (notice I'm not starting at steady state)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "Fin = 1\n",
    "h = 1\n",
    "V = A*h\n",
    "x0 = x = 0.7"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Valve characterisitic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(x):\n",
    "    return alpha**(x - 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Integration time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "ts = numpy.linspace(0, 100, 1000)\n",
    "dt = ts[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that I have reordered the equations here so that they can be evaluated in order to find the volume derivative."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "hs = []\n",
    "for t in ts:\n",
    "    h = V/A\n",
    "    Fout = K*f(x)*numpy.sqrt(h)\n",
    "    dVdt = Fin - Fout\n",
    "    V += dVdt*dt\n",
    "    \n",
    "    hs.append(h)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x10e93ca20>]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAGBRJREFUeJzt3X10XPV95/H3V6Nn2fKDJD9bCBvHgHkIjkiBdsHJpq0DNKQJ29TtNk03Od4kbNv0pG26p7thk5zTlj5tm5MmHC+lNDkJaRMIoXlomhBY04VAZDBgsA02fhJ+0Eiy9SyNRvPdP+bKlo2kka07vpp7P69z5szcub+Z+V6u/eHn3/3de83dERGReCmLugAREQmfwl1EJIYU7iIiMaRwFxGJIYW7iEgMKdxFRGJI4S4iEkMKdxGRGFK4i4jEUHlUP9zY2OgtLS1R/byISEnasWNHp7s3FWoXWbi3tLTQ1tYW1c+LiJQkMzs0k3YalhERiSGFu4hIDCncRURiSOEuIhJDCncRkRhSuIuIxJDCXUQkhiKb5y7J4O5kc87oWI7RMSc7lju9nB1zsrkcmWz+eeL6nDs5h5w77k4ux1nvjb/28de58bZn2o2d/qzjp+sJnifUd3a9E14Hrd78mbPXn/u5id977mcn+/xkbRIlgbf6bG1ZzM1vKXge0qwUDHczux+4Hehw96smWb8J+DZwIHjrYXf/bJhFysXj7vQMjXJycJTeoVF6h0fpHcoGz6P0Dedf949kGR4dY3g0x1BmjKHRsWA5/3ooM8ZwNkcmm4t6k6QEmEVdwcX10VvWRh/uwAPAF4AvT9PmSXe/PZSKpGgy2RzHeoY40j1E+8lB2k8O0dE3TGd/hnTfCJ39+cfo2NQ9qVSZUV9dTm1lOTWVKWoq8o/51eUsmV9FTWWK6vJU/rkiRWV5GRVlRkV5GeVlRkWqjPKUUVGWfy5PlVGZMsqD5YpUvl2qzDAzygzKzCgzw8Zfl515703ry868ZwYpO/M94wwbf5F/OnsRm5A0Z947+7OThdFUbWb0vUlLNym6guHu7tvNrKX4pUhYeodHee1EH6+e6OfVE328eqKPA+kBjvcOk5uQ26kyo6Gukqb5VTTOq2L9svk0zquicV4li+sqWVBTQX1NBfXVFdTXlFNfXUFtZUpBJFICwhpzv9HMXgCOAr/v7i+H9L1SQCabY/exXp4/fJKdR06x88gpDnYNnl5fU5Fi3dJ53LCmgVWLa1m9qIZVi2pZtaiG5QuqKU/pmLpIHIUR7s8Bl7h7v5ndCjwCrJusoZltBbYCNDc3h/DTyePu7E8P8ORraf79tU6efr2LwcwYAEvmV3Fd80L+U+tqLl82n7csnc/KhTWUlamnLZI0sw53d++d8Pp7ZvZFM2t0985J2m4DtgG0trYm7xD5BXJ3Xj7ay3dfOsZ3XzzG4e58z7yloZb3bVzJjWsaua55IcsXVGvIRESAEMLdzJYBJ9zdzezt5OfOd826MiHdN8I/tx3hG21HONg1SKrMuGltA1tvXsPN65pobqiNukQRmaNmMhXyQWAT0Ghm7cDdQAWAu98L3Al8zMyywBDwq37u5GE5L20Hu/mHpw7yg13HyeacG9Ys5qO3rOUXNixjcV1l1OWJSAmYyWyZLQXWf4H8VEmZBXfnqf1dfP6x13jmQDcLair40E0tbPmZZtY2zYu6PBEpMTpDdQ54qb2Hz33nFZ492M3S+io+ffuVbHl7MzWVqahLE5ESpXCPUEffMPd8fy8PPddOQ10ln7tjA79y/WqqyhXqIjI7CvcIuDvf3nmUux99maHMGP/1ljXc9Y7LqK+uiLo0EYkJhftF1tU/wqceeokf7T7BxuaF/Pmd13LZEo2pi0i4FO4X0Y5D3dz11efpHszwP267gt/62UtJ6QQjESkChftF8pWnD/KZf3mFlYtqePhjN3HVygVRlyQiMaZwL7JczvnT7+/m/zx5gHddsYS/+pW3sqBGY+siUlwK9yLKZHP83j/t5LsvHeNDN7XwP2+/UsMwInJRKNyLJJPN8fGvPsePdp/gj2+9go/8h0t13RcRuWgU7kWQyea462v5YP/sHRv44I0tUZckIgmji3mHLJdzPvmNF/jhKwp2EYmOwj1kf/Fve/mXF47yh5vXK9hFJDIK9xB97ZnDfOmJ/fzazzTzsVvWRl2OiCSYwj0kbQe7+fS3d7FpfROffc8GHTwVkUgp3EPQ2T/Cf/va86xcVMPnt1yn+5KKSOQ0W2aWcjnnE1/fycnBDA9//CZd/EtE5gSF+yzd//8O8O/7Ovmz913NhhW6pICIzA0aP5iFfR19/PkP9vKuK5bygetXR12OiMhpCvcLlB3L8cl/foG6yhR/8r6rdABVROYUDctcoC8/fYgX2nv4wq9dx5L51VGXIyJyFvXcL0BH7zB//cNX2bS+iduuXh51OSIib6JwvwB/8r3dZLI5/tcvaT67iMxNCvfz9OyBbh7ZeZSP3rKGlsa6qMsREZmUwv08uOdvvLGsvpqPbbos6nJERKakcD8P//bKCZ4/fIpPvGsdNZWpqMsREZmSwn2GsmM5/uIHe1nbVMedb1sVdTkiItNSuM/Qt55/g30d/fzBL67XtWNEZM5TSs3AWM754hP72bCinl/csCzqckREClK4z8C/7jrOgc4BPr7pMk19FJGSoHAvwN354hP7WNNYx+ar1GsXkdJQMNzN7H4z6zCzXQXaXW9mY2Z2Z3jlRW/7a528fLSXj96yllSZeu0iUhpm0nN/ANg8XQMzSwH3AD8IoaY55b4nX2dZfTXvvW5l1KWIiMxYwXB39+1Ad4Fmvw08BHSEUdRcsT/dz5OvdfKfb2imslwjWCJSOmadWGa2Evhl4N4ZtN1qZm1m1pZOp2f700X3lacPUZEyPnB9c9SliIiclzC6o38DfMrdxwo1dPdt7t7q7q1NTU0h/HTxDIxkeWhHO7ddvZym+VVRlyMicl7CuJ57K/D1YIpgI3CrmWXd/ZEQvjsy33r+DfpGsvzGjS1RlyIict5mHe7ufun4azN7APhOqQc7wNeeOcyGFfVsbF4YdSkiIuetYLib2YPAJqDRzNqBu4EKAHcvOM5einYf6+WVY7185j26XruIlKaC4e7uW2b6Ze7+oVlVM0c8tKOdipTxnmtXRF2KiMgF0fy+c4yO5Xhk5xu88/IlLKqrjLocEZELonA/x/ZX03T2Z3j/Rl3WV0RKl8L9HA89187iuko2rV8SdSkiIhdM4T5B/0iWH+3u4JeuWa4zUkWkpCnBJnh8TweZbI7brtGBVBEpbQr3Cb6/6xhN86t42yWLoi5FRGRWFO6BwUyWx/ek2bxhmS7tKyIlT+Ee+L970wyNjvHuq3VDDhEpfQr3wPd2HaehrpK3tyyOuhQRkVlTuAMj2TF+vPsEv7BhKeUp/ScRkdKnJAOePdDNQGaMn79yadSliIiEQuEO/HhPB1XlZdy4pjHqUkREQqFwB57Ym+bGtQ3UVKaiLkVEJBSJD/cDnQMc6BzgHbrcgIjESOLD/Ym9+Xt6K9xFJE4SH+6P702ztqmO5obaqEsREQlNosN9MJPlJ693qdcuIrGT6HD/6cGTZLI5bn5LU9SliIiEKtHh/tT+TipSxvU6K1VEYibR4f70/i6uW71IUyBFJHYSG+49Q6PseqOHG9c2RF2KiEjoEhvuzx7oJudwk8JdRGIoseH+1P5OqsrLeGvzwqhLEREJXWLD/en9XVzfspiqco23i0j8JDLcu/pH2HO8T+PtIhJbiQz3nx7sBuCGNQp3EYmnRIb7jkMnqSwv4+qVC6IuRUSkKBIb7tesXEBleSI3X0QSIHHpNjw6xq43ennbJYuiLkVEpGgSF+4vH+0hM5Zjo8JdRGKsYLib2f1m1mFmu6ZYf4eZvWhmO82szcx+Lvwyw7Pj0EkANjYr3EUkvmbSc38A2DzN+seAa939rcB/Ae4Loa6iee7QKZoX19I0vyrqUkREiqZguLv7dqB7mvX97u7BYh3gU7WNmruz4/BJjbeLSOyFMuZuZr9sZnuA75LvvU/VbmswdNOWTqfD+Onz0n5yiHTfiMbbRST2Qgl3d/+Wu18OvBf43DTttrl7q7u3NjVd/BtknBlv1/VkRCTeQp0tEwzhrDWzxjC/NywvtvdQXVHG+qXzoy5FRKSoZh3uZnaZmVnweiNQCXTN9nuLYdcbPVyxvJ7yVOJmgIpIwpQXamBmDwKbgEYzawfuBioA3P1e4P3AB81sFBgCPjDhAOuckcs5Lx/t4f1vWxV1KSIiRVcw3N19S4H19wD3hFZRkbzeOcBAZoyrdD0ZEUmAxIxP7HqjB0AXCxORREhMuL/Y3kNVeRnrlsyLuhQRkaJLTLjrYKqIJEkikm78YKqGZEQkKRIR7uMHUxXuIpIUiQj38YOpmikjIkmRiHB/+WgPleVlrFuqg6kikgyJCPc9x/t4y9J5VOhgqogkRCLSbs/xPtYvrY+6DBGRiyb24d7VP0K6b4QrlutiYSKSHLEP973H+wBYv0zhLiLJEftw36NwF5EEin247z3eR0NdJU3zdM9UEUmO2If7nhN9rF82n+CS8yIiiRDrcM/lnFeP92lIRkQSJ9bhfrh7kKHRMa5YpmmQIpIssQ53HUwVkaSKdbjv68iH+2W6hruIJEysw31/eoAVC6qpqyp4N0ERkViJdbjv6+hnrXrtIpJAsQ13d2d/up+1TQp3EUme2Ib7sZ5hBjNjGm8XkUSKbbjvT/cDqOcuIokU23Df15EPd/XcRSSJYh3u9dXlNM6rjLoUEZGLLrbhvj/dz2VL5umaMiKSSLEN930dAxpvF5HEimW49wyO0tk/ovF2EUmsWIb7/s78wdQ16rmLSELFMtwPdQ0AcGljXcSViIhEo2C4m9n9ZtZhZrumWP/rZvZi8HjKzK4Nv8zzc7BzEDNYvbgm6lJERCIxk577A8DmadYfAG5x92uAzwHbQqhrVg51DbBiQQ1V5amoSxERiUTByyW6+3Yza5lm/VMTFn8CrJp9WbNzoGtQQzIikmhhj7l/GPj+VCvNbKuZtZlZWzqdDvmnzzjUNcAlDbVF+34RkbkutHA3s3eQD/dPTdXG3be5e6u7tzY1NYX102c5NZjh1OAoLQ3quYtIcoVyFwszuwa4D3i3u3eF8Z0X6lDXIIB67iKSaLPuuZtZM/Aw8Bvu/ursS5qdg5oGKSJSuOduZg8Cm4BGM2sH7gYqANz9XuDTQAPwxeA6Lll3by1WwYWcmQapnruIJNdMZstsKbD+I8BHQqtolg51DbC8vprqCk2DFJHkit0Zqge7BmjRkIyIJFzswv1Q1yCXaKaMiCRcrMK9fyRL10BGM2VEJPFiFe5HuvPTIFcvUriLSLLFKtzbTw4BsGqRLhgmIskWq3A/3XPXNEgRSbhYhXv7ySFqK1Msqq2IuhQRkUjFLNwHWbWoRjfFFpHEi1W4Hzk5xCodTBURiVe4t58cZLUOpoqIxCfcewZH6RvOqucuIkKMwv3IyfxMGU2DFBGJUbiPz3HXNEgRkViFu3ruIiLjYhTuQ8yrKmdBjea4i4jEKNw1x11EZFxswv1I95CGZEREArEJ96Onhli5UOEuIgIxCfe+4VH6RrIsW6BwFxGBmIT78Z5hAFYsrI64EhGRuSEW4X40CPfl6rmLiAAxCfdjp/InMC1foJ67iAjEJdx7hjGDZQp3EREgNuE+RNO8KipSsdgcEZFZi0UaHusZ1pCMiMgEsQj3o6eGdDBVRGSCkg93d8/33DUNUkTktJIP997hLIOZMVao5y4iclrJh/uxnmAapHruIiKnFQx3M7vfzDrMbNcU6y83s6fNbMTMfj/8Eqd37JROYBIROddMeu4PAJunWd8N/A7wl2EUdL6O9ugEJhGRcxUMd3ffTj7Ap1rf4e4/BUbDLGymjp0apsxgyfyqKH5eRGROisGY+zBL66sp1wlMIiKnXdRENLOtZtZmZm3pdDqU7zzRmw93ERE546KGu7tvc/dWd29tamoK5Tvz4a4hGRGRiUp+LKOjb0Q9dxGRc5QXamBmDwKbgEYzawfuBioA3P1eM1sGtAH1QM7MPgFc6e69Ras6MDw6Rs/QqMJdROQcBcPd3bcUWH8cWBVaReeho3cEgCbNlBEROUtJD8t09OVPYFLPXUTkbCUd7ieCnrvmuIuInK2kw109dxGRyZV0uJ/oHaEiZSyqrYi6FBGROaWkw72jd5gl86sxs6hLERGZU0o73PtGWKITmERE3qSkw/1E77AOpoqITKKkw11np4qITK5kw11np4qITK1kwz3dp7NTRUSmUrLhfqJXc9xFRKZSwuGe77nrcr8iIm9WsuE+fnbqkvnquYuInKuEw11np4qITKVkw72zb4SGuiqdnSoiMonSDff+ERrnV0ZdhojInFSy4d41kKGhTgdTRUQmU7Lh3tk3QuM8hbuIyGRKMtzdnc7+jIZlRESmUJLh3jeSJTOWo1HDMiIikyrJcO8MLj2gnruIyORKM9z7MwAacxcRmUJJhntXf9BzV7iLiEyqJMO9Mwj3hnkalhERmUxJhnu6P4MZLK5VuIuITKYkw72rf4TFtZWUp0qyfBGRoivJdOzsH9GQjIjINEo03DM6mCoiMo2SDPeufl16QERkOiUZ7p39GQ3LiIhMo2C4m9n9ZtZhZrumWG9m9nkz22dmL5rZxvDLPGN4dIz+kax67iIi05hJz/0BYPM0698NrAseW4Evzb6sqaWDSw80KdxFRKZUMNzdfTvQPU2TO4Ave95PgIVmtjysAs/VNRBcekDXlRERmVIYY+4rgSMTltuD94pi/KJhulGHiMjUwgj3yW5i6pM2NNtqZm1m1pZOpy/oxxbWVrB5wzKWL6y+oM+LiCRBeQjf0Q6snrC8Cjg6WUN33wZsA2htbZ30fwCFtLYsprVl8YV8VEQkMcLouT8KfDCYNXMD0OPux0L4XhERuUAFe+5m9iCwCWg0s3bgbqACwN3vBb4H3ArsAwaB3ypWsSIiMjMFw93dtxRY78BdoVUkIiKzVpJnqIqIyPQU7iIiMaRwFxGJIYW7iEgMKdxFRGLI8pNdIvhhszRw6AI/3gh0hlhOKdA2J4O2ORlms82XuHtToUaRhftsmFmbu7dGXcfFpG1OBm1zMlyMbdawjIhIDCncRURiqFTDfVvUBURA25wM2uZkKPo2l+SYu4iITK9Ue+4iIjKNkgt3M9tsZnuDG3L/UdT1FIOZrTazx81st5m9bGa/G7y/2Mx+aGavBc+Loq41TGaWMrPnzew7wfKlZvZMsL3/ZGaxureimS00s2+a2Z5gX9+YgH38e8Gf6V1m9qCZVcdtP5vZ/WbWYWa7Jrw36X4NLpX++SDPXjSzjWHVUVLhbmYp4O/I35T7SmCLmV0ZbVVFkQU+6e5XADcAdwXb+UfAY+6+DngsWI6T3wV2T1i+B/jfwfaeBD4cSVXF87fAv7r75cC15Lc9tvvYzFYCvwO0uvtVQAr4VeK3nx8ANp/z3lT79d3AuuCxFfhSWEWUVLgDbwf2ufvr7p4Bvk7+Bt2x4u7H3P254HUf+b/0K8lv6z8Gzf4ReG80FYbPzFYBtwH3BcsGvBP4ZtAkbttbD9wM/D2Au2fc/RQx3seBcqDGzMqBWuAYMdvP7r4d6D7n7an26x3Alz3vJ8BCM1seRh2lFu4X9Wbcc4GZtQDXAc8AS8fvchU8L4mustD9DfCHQC5YbgBOuXs2WI7bvl4DpIF/CIai7jOzOmK8j939DeAvgcPkQ70H2EG89/O4qfZr0TKt1MJ9xjfjjgMzmwc8BHzC3XujrqdYzOx2oMPdd0x8e5KmcdrX5cBG4Evufh0wQIyGYCYTjDPfAVwKrADqyA9LnCtO+7mQov05L7Vwn/HNuEudmVWQD/avuvvDwdsnxv/JFjx3RFVfyH4WeI+ZHSQ/1PZO8j35hcE/3yF++7odaHf3Z4Llb5IP+7juY4B3AQfcPe3uo8DDwE3Eez+Pm2q/Fi3TSi3cfwqsC46uV5I/GPNoxDWFLhhv/ntgt7v/9YRVjwK/Gbz+TeDbF7u2YnD3/+7uq9y9hfw+/bG7/zrwOHBn0Cw22wvg7seBI2a2PnjrPwKvENN9HDgM3GBmtcGf8fFtju1+nmCq/foo8MFg1swNQM/48M2suXtJPcjfjPtVYD/wx1HXU6Rt/Dny/zR7EdgZPG4lPw79GPBa8Lw46lqLsO2bgO8Er9cAz5K/+fo3gKqo6wt5W98KtAX7+RFgUdz3MfAZYA+wC/gKUBW3/Qw8SP6Ywij5nvmHp9qv5Idl/i7Is5fIzyQKpQ6doSoiEkOlNiwjIiIzoHAXEYkhhbuISAwp3EVEYkjhLiISQwp3EZEYUriLiMSQwl1EJIb+P0I8V1a6yBrWAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(ts, hs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PI Control\n",
    "\n",
    "\n",
    "Now we can include a controller controlling the level by manipulating the valve fraction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.signal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's do a PI controller:\n",
    "\n",
    "$$ G_c = K_c(1 + \\frac{1}{\\tau_I s}) = \\frac{K_C \\tau_I s + K_Cs^0}{\\tau_I s + 0s^0} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "Kc = -1\n",
    "tau_i = 5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In versions of scipy < 1.0, `Gc` would automatically have a `.A` attribute. After 1.0, we need to convert to state space explicitly with `.to_ss()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "Gc = scipy.signal.lti([Kc*tau_i, Kc], [tau_i, 0]).to_ss()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "hsp = 1.3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "V = 1\n",
    "hs = []\n",
    "xc = numpy.zeros([Gc.A.shape[0], 1])\n",
    "for t in ts:\n",
    "    h = V/A\n",
    "    e = hsp - h  # we cheat a little here - the level we use to calculate the error is from the previous time step\n",
    "    \n",
    "    # e is in the input to the controller, yc is the output\n",
    "    dxcdt = Gc.A.dot(xc) + Gc.B.dot(e)\n",
    "    yc = Gc.C.dot(xc) + Gc.D.dot(e)\n",
    "    \n",
    "    x = x0 + yc[0,0]  # x0 is the controller bias\n",
    "    \n",
    "    Fout = K*f(x)*numpy.sqrt(h)\n",
    "    dVdt = Fin - Fout\n",
    "    \n",
    "    V += dVdt*dt\n",
    "    xc += dxcdt*dt\n",
    "    \n",
    "    hs.append(h)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1c1a0f31d0>]"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAFtRJREFUeJzt3X2QXXWd5/H3t7vTnUeTkHQACSE8hPiAMEJEQUVQdwRmVrTKWWWsHZ3FSc2u46hDuaO1W7rr7tbUPOzu6Iq6DCLq7oCIllIUyloMM6AjSFBgkMcMCESQBEgCSbpv33v7u3/c20k/3NvdSe5Nc07er6que8+5J/d+T53kk19/z++eE5mJJKlceua6AElS5xnuklRChrsklZDhLkklZLhLUgkZ7pJUQoa7JJWQ4S5JJWS4S1IJ9c3VB69cuTLXrl07Vx8vSYV01113PZuZgzNtN2fhvnbtWjZt2jRXHy9JhRQRj89mO9syklRChrsklZDhLkklZLhLUgkZ7pJUQoa7JJWQ4S5JJVTqcN/2YoWv/Ogx/uHhbXg7QUmHkzn7ElO3vThc5aIv/Iindg4D8IGzjuMz//LV9PTEHFcmSd1X2nD/m9se46mdw3z935zJrQ9v44ofPcbyRf187O0nz3VpktR1pQz3zOR7d/+KN69byTknD/LmdSt5fvcIn7v5EU47dhnnrV811yVKUleVsuf+8DO7ePy5PVxwytEARAT/7d2vYf2RS7j02nv4dbNVI0llNWO4R8SVEbE1Iu6bYbvXRUQ9It7TufIOzM+f2A7AWSeu2LtuQX8vX/jd0xkaqfOxb/6c+qgnWCWV12xG7lcB50+3QUT0An8O3NSBmg7aPVt2sHTBPNauWDhh/UmrFvPZi17N7Y8+zxf+bvMcVSdJ3TdjuGfmrcDzM2z2EeDbwNZOFHWw/ulXOzl19VIips6Mec8Zq3nXb7ycz938MHc8+twcVCdJ3XfQPfeIOAZ4N/Dlgy/n4GUmj27bzUmrFrd8PSL4r+9+DWuOWMhHr7mb53ePHOIKJan7OnFC9a+BP83M+kwbRsTGiNgUEZu2bdvWgY+eauuLFfaM1Dl+5aK22ywe6OMLv3s6z+8e4RPfuscvOEkqnU6E+wbgmoj4JfAe4IsR8a5WG2bm5Zm5ITM3DA7OeJeoA/LYs7sBWLuifbgDnHLMUj514Su4+cGtfPXHv+xKLZI0Vw56nntmHj/2PCKuAm7IzO8e7PseqF82w326kfuYD569lh9vfo4/+/4DnHbsUs447ohulydJh8RspkJeDfwEWB8RWyLikoj4w4j4w+6Xt/+ebs5hP2rp/Bm3jQj+6ndOZfXyhVzytU1s3rqr2+VJ0iExm9kyF2fm0Zk5LzNXZ+ZXMvPLmTnlBGpmfjAzr+tOqbOz9cUKKxf3M693dh2nZQv7+drvn0lfT/CBK3/KMy/4BSdJxVe6b6hufWGYwSUzj9rHW7NiIV/94Jns2DPCe//3T3hqx1CXqpOkQ6N04f7Mi8Mc+bKB/f5zr1m9lK9f8nqe2zXCey//CU8+v6cL1UnSoVG6cN/6QoUj93PkPuaM45bzfz70el4YqvHuL/6Yux6f6btbkvTSVKpwr48mz+6qHNDIfcxpxy7j2//2bBYN9HHx5XfwnZ9t6WCFknRolCrcn9tVYTRh8GUHNnIfc9KqxXz3372R049bxp9cew+f+NY97KrUOlSlJHVfqcL9mRcqABy55MBH7mOWL+rnG5e8no+89SS+/bMt/Nbnb/NaNJIKo2Th3pjGuOogR+5j5vX2cOlvrueajWdRH03ee/ntfPybd7P1RadLSnppK1W4b9/TuAjYikX9HX3fM48/gh9+/C18+LwTueHepzj3L/+ev/jBg2z3omOSXqJKFe47h6oALFs4r+PvvaC/l0+84xXc9LFzeNsrj+RL//DPvPkvbuHPvv+A0yYlveSUKtx37KnS2xMsHujerWFPGFzM/7r4tfzgo+fwlpMHueK2xzjnL2/hkqvu5If3P0OlNuPFMSWp60p1g+wdQyMsWzCv5U06Om39UUu47P2n8/TOIf72jie4+qdPcPODW1kyv493vPoo3vHqozjrxBVd/Y9GktopVfLs2FNl6YLOt2Smc/TSBVz6m+v547et40ebn+WGe57mpvt+zXV3baGvJzh9zXLevG4lZ6xdzqmrlxn2kg6JUiXNzqEqS7vQb5+Neb09nLd+FeetX0Wldgp3Pb6d2x55ltse2cZ//+HDAETAulWLOXX1MtYfuYSTVi3mpFWLOWbZAnp6uv/bhqTDR6nCfceeKisXd3amzIEY6Ovl7BNXcvaJK/nT81/B9t0j3LNlB3c/uYN7ntzB3z+0levu2jJu+x7WrljEy5fN5+hlCzhm2QKOXjqfo5bOZ8WiAZYvmsfyhbO/0qUklSvch0ba3jt1Li1f1M+561dx7vpVe9dt3z3C5m272Ly18fP4c3t4eucQdz+5g+17qi3fZ8n8Po5Y1M/yhf0smd/Hov4+Fg30sWigt/HYP/bYx8C8Hvp7e5jX20N/37if3omPfb1BbwS9PUFPz7jnex85JOcwJHVWucJ9DnruB2r5on5et+gIXrd26t2fhkbqPL1ziF+/MMz23VWe3zPC9t0jPL97hO17Go+7KjWeeWGY3ZU6u0dq7K7UqNa7cy/YCOiNyeHP3uVG9sfebWPcn4uW62PCe48tBjHu+b7txr39lPVl/n8nKOfOlfmYzdbvbDiWS950/MwbHoTShHutPsqLw7WuzHE/1Bb093LC4GJOGNy/30JGaqPsrtTYVakxUh9lpDZKtfk4UhulUh+lWhvd+9rY66PZuOjaaCb10aSeyehotl1fH4XR3Ld+7L+Uxn3Gc+/zsfuOJznu+cT1TFifLbaZuJ4W71lGZd23fX9bymV/j9fyQ5BTpQn3F4YbF/ZaVpCRezc0Wi/9LO/wN3QlFU9pztDtaF56YNlCg02SyhPuzUsPFKXnLkndVJpwf7HZlnnZgtJ0miTpgJUm3Hc1w33xgCN3SSpPuFcabZlFA71zXIkkzb3ShPtYW2aJI3dJKk+47640LrXryF2SShTuuypV5s/roc/rr0hSmcK95slUSWoqUbjXWTLfaZCSBGUK9+GqN8KQpKbyhHul5slUSWoqUbjX7blLUlOJwr1qz12SmsoT7sO2ZSRpzIzhHhFXRsTWiLivzevvj4h7mz//GBGndb7MmTkVUpL2mc3I/Srg/Glefwx4S2aeCvwX4PIO1LVfKrU61XralpGkphnTMDNvjYi107z+j+MWbwdWH3xZ+2ffFSENd0mCzvfcLwG+3+H3nNGuSiPcFxnukgR08B6qEXEejXB/0zTbbAQ2AqxZs6ZTH73vomH9nlCVJOjQyD0iTgWuAC7KzOfabZeZl2fmhszcMDg42ImPBmCo2hi5LzDcJQnoQLhHxBrgO8C/zsyHD76k/bdnpDFyX9hvW0aSYBZtmYi4GjgXWBkRW4DPAPMAMvPLwKeBFcAXIwKglpkbulVwK/vC3ZG7JMHsZstcPMPrHwI+1LGKDsBQM9xty0hSQym+oerIXZImKkm4N06oLpxnz12SoCThbltGkiYqRbjvqdbp6wn6+0qxO5J00EqRhkMjdUftkjROKcJ9z0jNk6mSNE5Jwr3uF5gkaZxShPvQSJ0F8xy5S9KYUoR7Y+RuuEvSmHKEe9UTqpI0XinCfcgTqpI0QSnC3ROqkjRRKcJ92LaMJE1QinDfM1JnobNlJGmvwod7ZjJUdbaMJI1X+HAfro6SCQvsuUvSXoUP972X+3XkLkl7lSDcvdyvJE1W+HAfqnoXJkmarPDh7i32JGmqEoR7o+e+wFvsSdJehQ/3IUfukjRF4cPdtowkTVX4cB87oTrfb6hK0l6FD/fhqlMhJWmy0oS7I3dJ2qcE4T4KwPy+wu+KJHVM4RNxuFqnryfo6y38rkhSxxQ+EYeqdVsykjRJ4cN9uDrK/HmF3w1J6qjCp2LFkbskTVH4cB+uGe6SNFnxw922jCRNUfhUHBqpM7/PkbskjTdjuEfElRGxNSLua/N6RMTnI2JzRNwbEad3vsz2bMtI0lSzGblfBZw/zesXAOuaPxuBLx18WbPXaMsY7pI03ozhnpm3As9Ps8lFwNez4XZgWUQc3akCZ9KYLVP47pIkdVQnUvEY4Mlxy1ua66aIiI0RsSkiNm3btq0DH934hqojd0maqBPhHi3WZasNM/PyzNyQmRsGBwc78NEwXHO2jCRN1olU3AIcO255NfBUB953VpwtI0lTdSLcrwd+rzlr5g3Azsx8ugPvO6PMZLhW91rukjTJjHeVjoirgXOBlRGxBfgMMA8gM78M3AhcCGwG9gC/361iJxupj5LptdwlabIZwz0zL57h9QQ+3LGK9sPYtdwHvJa7JE1Q6FSseBcmSWqp0OHuzbElqbVCh/veW+w5FVKSJih0Ko7dHHuBI3dJmqAU4W5bRpImKna412zLSFIrhU7FoZHGyH3Ab6hK0gSFDvdKzbaMJLVS6HDfe0LVyw9I0gQFD/dmz91vqErSBIVORWfLSFJrBQ/3sdkyhrskjVfocB+q1pnXG/T2tLpfiCQdvgod7sNVb9QhSa0UOtwrtTrznSkjSVMUOtyHq94/VZJaKXQy2paRpNYKHe5D1bozZSSphUKH+3C1bltGkloodDI2eu6O3CVpsoKHu20ZSWql0OFeqTlyl6RWCh3ujdkyhd4FSeqKQiejI3dJaq3Q4T5crTPgyF2Spih0MlZqoww4FVKSpihsMtbqo9RH0/unSlILhQ33Sq1xLXfbMpI0VWGT0XCXpPYKm4yVWuMWewPOlpGkKYob7ntvsVfYXZCkrilsMu5ryzhyl6TJChzuzbaMPXdJmmJWyRgR50fEQxGxOSI+2eL1NRFxS0T8PCLujYgLO1/qRI7cJam9GcM9InqBy4ALgFcBF0fEqyZt9h+BazPztcD7gC92utDJxnrufolJkqaaTTKeCWzOzEczcwS4Brho0jYJvKz5fCnwVOdKbG24altGktrpm8U2xwBPjlveArx+0jb/Cfh/EfERYBHw9o5UNw3bMpLU3myGvdFiXU5avhi4KjNXAxcC34iIKe8dERsjYlNEbNq2bdv+VzuOJ1Qlqb3ZJOMW4Nhxy6uZ2na5BLgWIDN/AswHVk5+o8y8PDM3ZOaGwcHBA6u4ae/I3Z67JE0xm2S8E1gXEcdHRD+NE6bXT9rmCeBtABHxShrhfnBD8xlU9vbcbctI0mQzhntm1oA/Am4CHqAxK+YXEfHZiHhnc7NLgT+IiHuAq4EPZubk1k1HeW0ZSWpvNidUycwbgRsnrfv0uOf3A2/sbGnTM9wlqb3CJmOlVqevJ+jrLewuSFLXFDYZK9VRR+2S1EZh07Fxiz1PpkpSKwUOd2+OLUntFDYdh23LSFJbhU3HxsjdtowktVLgcB/126mS1EZh09HZMpLUXmHT0baMJLVX4HB35C5J7RQ2HSu1UeY7z12SWipwuDvPXZLaKWw6VqrOlpGkdgqbjo2eu20ZSWqlwOFuW0aS2ilkOmamlx+QpGkUMh1H6mP3T7UtI0mtFDLcvQuTJE2vkOlYqRrukjSdQqZjpVYHcLaMJLVR0HAf67kXsnxJ6rpCpuO+towjd0lqpZjhPtaWceQuSS0VMh2dLSNJ0ytkOu4Ld9syktRKMcO9OjZbppDlS1LXFTIdh5sj9/n23CWppUKm476Ru20ZSWqlmOHuCVVJmlYh09ETqpI0vYKGu/PcJWk6hUzHsW+o9vcWsnxJ6rpCpmOlNkp/Xw89PTHXpUjSS9Kswj0izo+IhyJic0R8ss02/yoi7o+IX0TE33a2zIm8xZ4kTa9vpg0iohe4DPgXwBbgzoi4PjPvH7fNOuBTwBszc3tErOpWweDNsSVpJrMZ/p4JbM7MRzNzBLgGuGjSNn8AXJaZ2wEyc2tny5yo4v1TJWlas0nIY4Anxy1vaa4b72Tg5Ij4cUTcHhHnd6rAViq1ujNlJGkaM7ZlgFZnLbPF+6wDzgVWA7dFxCmZuWPCG0VsBDYCrFmzZr+LHTNctS0jSdOZzfB3C3DsuOXVwFMttvleZlYz8zHgIRphP0FmXp6ZGzJzw+Dg4IHW7AlVSZrBbBLyTmBdRBwfEf3A+4DrJ23zXeA8gIhYSaNN82gnCx2vcULVcJekdmZMyMysAX8E3AQ8AFybmb+IiM9GxDubm90EPBcR9wO3AJ/IzOe6VXSlNsrAPNsyktTObHruZOaNwI2T1n163PME/qT503WVap2BJQOH4qMkqZAK2dsYsS0jSdMqZEJWaqPMty0jSW0VNNydLSNJ0ylkQlac5y5J0ypmuNdG/YaqJE2jcAk5OpqM1D2hKknTKVxCeos9SZpZAcO9eYs9R+6S1FbhEnLvyN2euyS1VbiEHLt/qm0ZSWqveOFuW0aSZlS4hNx3QrVwpUvSIVO4hBwbuXv5AUlqr3jhXnXkLkkzKVxC7pst48hdktopYLh7QlWSZlK4hBxcMsCFrzmKZQvnzXUpkvSSNas7Mb2UnHHcEZxx3BFzXYYkvaQVbuQuSZqZ4S5JJWS4S1IJGe6SVEKGuySVkOEuSSVkuEtSCRnuklRCkZlz88ER24DHD/CPrwSe7WA5ReA+Hx7c58PDwezzcZk5ONNGcxbuByMiNmXmhrmu41Bynw8P7vPh4VDss20ZSSohw12SSqio4X75XBcwB9znw4P7fHjo+j4XsucuSZpeUUfukqRpFC7cI+L8iHgoIjZHxCfnup5uiIhjI+KWiHggIn4RER9trj8iIn4YEY80H5fPda2dFBG9EfHziLihuXx8RNzR3N9vRkT/XNfYSRGxLCKui4gHm8f6rMPgGH+8+Xf6voi4OiLml+04R8SVEbE1Iu4bt67lcY2Gzzfz7N6IOL1TdRQq3COiF7gMuAB4FXBxRLxqbqvqihpwaWa+EngD8OHmfn4SuDkz1wE3N5fL5KPAA+OW/xz4n8393Q5cMidVdc/ngB9k5iuA02jse2mPcUQcA/wxsCEzTwF6gfdRvuN8FXD+pHXtjusFwLrmz0bgS50qolDhDpwJbM7MRzNzBLgGuGiOa+q4zHw6M3/WfP4ijX/0x9DY1681N/sa8K65qbDzImI18FvAFc3lAN4KXNfcpGz7+zLgHOArAJk5kpk7KPExbuoDFkREH7AQeJqSHefMvBV4ftLqdsf1IuDr2XA7sCwiju5EHUUL92OAJ8ctb2muK62IWAu8FrgDODIzn4bGfwDAqrmrrOP+Gvj3wGhzeQWwIzNrzeWyHesTgG3AV5utqCsiYhElPsaZ+Svgr4AnaIT6TuAuyn2cx7Q7rl3LtKKFe7RYV9rpPhGxGPg28LHMfGGu6+mWiPhtYGtm3jV+dYtNy3Ss+4DTgS9l5muB3ZSoBdNKs898EXA88HJgEY22xGRlOs4z6drf86KF+xbg2HHLq4Gn5qiWroqIeTSC/f9m5neaq58Z+5Wt+bh1rurrsDcC74yIX9Jotb2Vxkh+WfPXdyjfsd4CbMnMO5rL19EI+7IeY4C3A49l5rbMrALfAc6m3Md5TLvj2rVMK1q43wmsa55d76dxMub6Oa6p45r95q8AD2Tm/xj30vXAB5rPPwB871DX1g2Z+anMXJ2Za2kc07/LzPcDtwDvaW5Wmv0FyMxfA09GxPrmqrcB91PSY9z0BPCGiFjY/Ds+ts+lPc7jtDuu1wO/15w18wZg51j75qBlZqF+gAuBh4F/Bv7DXNfTpX18E41fze4F7m7+XEijD30z8Ejz8Yi5rrUL+34ucEPz+QnAT4HNwLeAgbmur8P7+hvApuZx/i6wvOzHGPjPwIPAfcA3gIGyHWfgahrnFKo0RuaXtDuuNNoylzXz7J9ozCTqSB1+Q1WSSqhobRlJ0iwY7pJUQoa7JJWQ4S5JJWS4S1IJGe6SVEKGuySVkOEuSSX0/wEuARU9TEWsdQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(ts, hs)"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}